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1.  INTRODUCTION 


There  are  any  number  of  situations  where  one  would  like  to  obtain 
information  about  the  behavior  of  the  crosswind  along  a  path.  One  method 
of  remote  sensing  that  offers  promise  is  that  of  optical  scintillations 
using  correlation  techniques. 

Provided  that  the  medium  is  weakly  turbulent,  then  it  is  possible 
to  write  down  a  nonlinear  integral  equation  between  the  unknown  cross- 
wind  profile  and  the  space-time  covariance  function  of  the  log-amplitude 

of  the  incident  laser  radiation.^  Lee  and  Uarp^  have  shown  how  to  con¬ 
vert  the  nonlinear  integral  equation  into  a  linear  integral  equation, 
and  working  only  with  the  slope  at  zero  time  delay  of  the  space-time 
covariance  function  of  the  signal. 

Lawrence,  Ochs,  and  Clifford,^  in  an  important  paper,  developed  an 
experimental  procedure,  utilizing  the  linear  integral  equation  approach, 
whereby  they  measured  an  averaged  crosswind,  the  average  taken  with  re¬ 
spect  to  the  path-weighting  function.  Furthermore,  they  are  able  to  make 
measurements  effectively  in  real  time.  In  view  of  their  success,  it  is 
now  an  opportune  time  to  consider  the  inversion  problem  of  reconstruct¬ 
ing  the  crosswind  profile  itself  via  measurements  of  the  correlation 
slope. 

The  solution  of  the  linear  integral  equation,  Eq.  (2)  or  Eq.  (9), 
is  an  inverse  problem.  Inverse  problems  are  known  to  be  ill-posed  (i.e., 
numerically  unstable)  and  any  inversion  method  must  be  capable  of  a  rea- 

3  4 

sonably  robust  inversion  in  the  presence  of  noise.' ’  We  propose  to  use 
the  method  of  singular  value  decomposition  to  invert  the  integral  equa¬ 
tion  in  the  presence  of  noisy  input  data. 

There  have  been  previous  inversion  studies  of  Eq .  (2).  Peskoff^ 
carried  out  an  analytical  inversion  of  the  linear  integral  equation 
which,  although  mathematically  correct,  suffers  from  the  fact  that  his 
path  must  be  infinite  in  order  to  perform  certain  of  the  analytical 
manipulations.  A  second  attempt  to  invert  the  integral  equation  was 

^R.W.  Lee  and  J.C.  Harp,  "Weak  Scattering  in  Random  Media,  With  Applica¬ 
tions  to  Remote  Probing, "  Proa.  IEEE,  57,  375-406. 

2 

R.S.  Lawrence,  G.R.  Ochs,  and  S.F.  Clifford,  "Use  of  Scintillations  to 
Measure  Average  Wind  Across  a  Light  Beam,"  Appl.  Opt.,  11,  1972,  239-243. 

'  M.M.  Laventiev,  Some  Improperly  Posed.  Problems  of  Mathematical  Physics, 
(Springer-Verlag,  New  York,  1967). 

4 

A.N.  Ttkhonov  and  V.Y.  Aresenin,  Solutions  of  Ill-Posed  Problems, 
(Halsted  Press,  New  York,  1977). 

5 A.Peskoff ,  "Theory  for  Remote  Sensing  of  Wind-Velocity  Profiles ,"  Proc. 
IEEE,  59,  1971,  324-325. 


made  by  Shen,  using  some  naive  numerical  techniques.  Unfortunately, 
he  did  not  realize  that  the  problem  was  ill-posed;  consequently,  his 
results  and  conclusions  are  misleading.  A  substantial  advance  was  made 

n 

by  Heneghan  and  Ishimaru'  who,  recognizing  that  the  problem  was  ill- 
posed,  employed  an  inversion  scheme  dependent  on  statistical  regulariza- 

g 

tion.  The  relation  between  their  method/results  and  our  method/results 
is  briefly  discussed  in  Section  IV. 


II.  PRELIMINARIES 


The  space  time  covariance  function  and  the  crosswind  profile  are 
related  to  each  other  by  a  complicated  nonlinear  integral  equation  whose 
explicit  form  we  need  not  quote,  see  Eq.  5  of  Ref.  2.  Following  the 

suggestion  in  Ref.  1,  the  nonlinear  integral  equation  can  be  made  linear 
(.in  the  velocity  profile)  by  differentiating  it  with  respect  to  the  time 
covariance  function  and  then  setting  the  time  delay  equal  to  zero. 


The  resulting  linear  integral  equation  relating  the  slope  of  the 
space-time  covariance  function  at  zero  time  lag  and  the  crosswind  com¬ 


ponent  as  a  function  of  the  separation  between  detectors  is‘ 
'L 


E(p)  = 


/ 


0  dzC“(z)W(Z,p)  v(z) 

/ 1  dzC“(z)[z(l-z)]S/(> 
/n  n 


(1 


where 


v(z) 


distance  along  sight  path  (0 , L ) 
crosswind  component 

refractive  index  structure  coefficient 


p  =  distance  between  detectors 

H(.P)  =  slope  with  respect  to  time  of  space  time  covariance  func¬ 
tion  at  zero  time  lag 
W(z,p)=  path-weighting  function. 

See  Fig.  1  for  schematic  of  the  geometry.  We  caution  the  reader  that 
this  expression  is  valid  only  in  a  weakly  turbulent  medium. 


L.C.  Sheri,  "Remote  Probing  of  Atmospheric  and  Wind  Velocity  By  Milli¬ 
meter  Waves,"  IEEE  Trans.  Antennas  and.  Prop.,  AP-18,  1970,  493-497. 

7 

J.M.  Heneghan  and  A.  Ishimaru,  " Remote  Determination  of  the  Profiles  of 
the  Atmospheric  Structure  Constant  and  Wind  Velocity  Along  a  Line  of 
Sight  Path  by  a  Statistical  Inversion  Procedure , "  IEEE  Trans.  Antennas 
and  Prop.,  AP-22,  1974,  157-464. 

Q 

4.N.  Franklin,  "Well-Posed  Stochastic  Extensions  of  Ill-Posed  Linear 
Problems,"  J.  Math.  Anal.  Appl. ,  31,  1970,  682-716. 
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We  are  particularly  interested  in  the  path  uniform  case  [C^(z)  = 

constant]  characteristic  of  horizontal  propagation.*  Equation  (1J  then 
becomes 


Jo 


dzW(z,p)v(z) 


E(P)  = 


/ 


L  dz[z(l-z)]5/6 
0 


(2) 


The  path-weighting  function  for  a  spherical  wave,  assuming  that 
the  refractive  index  size  spectrum  uses  the  inertial  subrange  assump- 
9 

tion,  takes  the  form 


W(z.p)  '  (2.33)(kL)5/6  fQ  dK  K*573  sin“ 

' 2J j  (ZDK/2L) 
x  (ZDK/2L) 

where  k  =  wavenumber  of  laser  radiation  and  I)  is  the  diameter  of  the  two 
detectors.  The  basic  point  to  emphasize  is  that  the  refractivity  spec¬ 
trum  is  taken  as  proport ional  to  K" ^  ^ .  The  degree  to  which  departures 
from  the  power  11/3  influence  the  final  results  is  not  undertaken  in 
this  paper,  but  it  should  be  kept  in  mind  that  11/3  is  a  useful  working 
number  and  not  an  axiomatic  statement. 


At  this  stage  it  is  convenient  to  convert  to  dimensionless  vari¬ 
ables: 


s  =  K(AL 


1/2 


=  KF 


6  = 


a  - 


(AL) 

D 


1/2 


(AL) 


1/2 


P_ 

F 

o 

F 


(4) 


Here 

F  S  (AL) 1/2  (5) 

is  the  Fresnel  zone  number  which  has  the  dimensions  of  a  length. 

*See,  liowever,  the  qualifying  remarks  in  Sea.  IV. 

9 

V.I.  Tatarskii ,  "The  Effects  of  the  Turbulent  Atmosphere  on  Wave  Propa¬ 
gation (US  Dept,  of  Commerce,  Springfield,  VA ,  1971). 
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The  path-weight  i njj  function  W(z,p)  can  be  rewritten  as 


c/h  ..7/3 

WIZ.B)  =  (2.33)  (2TT):,/t>  (2.6^ 


(b) 


where  W^(Z,6)  is  dimensionless 


i(Z-8)  ■  f’  dss'5/3  si,,J  p-ST1  [■(*:./.■)  ’ 


(82s)  .  (7) 


The  denominator  of  Eq.  2  can  be  evaluated  in  terms  of  gamma  func 


t  ions 


:[z(l-z)]5/6  -  L8/3 

'o  r  f1-! 


t  dz| 

*n 


Thus  the  dimensionless  form  of  Eq.  12)  is 
E(B)  =  £  f0  dzK1(Z.8)v(Z) 


(8) 


(9) 


=  48.8938 


(10) 


The  integral  on  the  right-hanu  size  has  dimensions  (length/time)  since 

v(Z)  is  a  velocity.  Hence,  E^B)  has  dimension  (time)  *  as  it  should. 

W'e  caution  the  reader  that  E(B)  and  v(Z)  must  be  measured  in  the  same 
time  units.  Equation  (9)  is  an  integral  equation  of  the  first  kind  for 
the  unknown  v(Z)  in  terms  of  the  measured  E(B)  and  known  W(Z,8). 


III.  INVERSION  VIA  SINGULAR  VALUES 

Given  this  basic  preliminary  information,  we  now  pass  on  to  the  in¬ 
version  of  Eq.  (9)  using  the  method  of  singular  value  decomposition. 

The  integral  equation,  Eq.  (9),  was  first  discretized  using  N  point 
Simpson's  rule 

J.  Hn"(8m'Zn)v<Z„>  '  E<8m>  "11 

n=  l 


11 


"here  Hn  are  the  weight  factors  ami  the  quadrature  points.  Since  we 

generally  want  more  measured  values  than  reconstruction  points,  N,  we 

let  m  *  l,  2 . M  where  M  _>  N.  Mien  tiq.  ill)  is  converted  to  matrix 

from,  we  have 


where  II  is  an  M  x  M  orthogonal  matrix  (IHl"  =  IMl  =  1  ),  \  is  an  \  x  \ 

orthogonal  matrix.  Hie  matrix  '  is  M  \  \  with  non-negative  elements  on 
the  main  diagonal  and  zeros  elsewhere 


N 


I  0  | 

I  he  's  are  t  >rmed  the  singular  values  of  \  and  are  the  non-negative 

square  roots  i  the  eigenvalues  of  A  \  (this  is  the  mathematical  defini¬ 
tion  of  the  i'«,  but  they  are  never  evaluated  from  the  definition). 

12 


The  solution  becomes  clearer  if  the  right-hand  side  of  Eq.  (15)  is 
written  out  explicitly. 

^C.B.  Moler  and  G.E.  Forsyth,  Computer  Solution  of  Linear  and  Alge¬ 
braic:  Systems,  (Prentioe-Hall,  Engleuood  Cliffs,  NJ ,  1967). 

H C.L .  Lawson  and  R.J.  Hanson,  Solving  Least  Squares  Problems , 
(Prentice- Hall,  Englewood  Cliffs,  NJ,  1974). 
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k  <  N 


(19) 


u  +  b  „ 
n 


where  u  and  v  denote  the  n  column  vectors  of  U  and  V  corresponding 
n  n 

to  o  .  The  sum  only  runs  over  n  =  1  to  k,  where  o.  is  the  smallest 
n  k 

nonzero  singular  value.  Equation  (19)  shows  that  the  matrix  x  of  rank 

k  (j-N)  is  a  linear  combination  of  k  matrices  of  rank  one. 


The  ill-posed  nature  of  the  inversion  is  directly  evident.  The 
smaller  singular  values  entering  into  the  denominator  tend  to  greatly 
magnify  any  error  in  the  measured  data  vector  b,  thereby  resulting  in  a 
spurious  solution.  To  alleviate  tins,  the  expansion  is  terminated  be¬ 
fore  the  contamination  due  to  the  numerically  small  singular  values  sets 

1 

in.  One  way  to  achieve  this,  Blake  and  Barakat,  “  is  to  to  set 


where  the  criterion  for  picking  e  is 


— -  >>  noise.  (21) 

°0 

Physically  this  procedure  has  the  effect  of  ignoring  high  frequency  com¬ 
ponents  of  b  which  are  the  main  cause  of  the  numerical  instability. 
Nevertheless  this  procedure  entails  a  somewhat  arbitrary  judgment  as  to 
when  to  terminate  the  summation.  Unfortunately  so  do  any  other  criteria. 

An  added  advantage  to  the  use  of  singular  values  is  that  the  singu¬ 
lar  values  are  very  stable  to  perturbations  in  the  matrix  elements,  in 
that  perturbations  of  the  matrix  elements  produce  perturbations  in  the 
singular  values  of  the  same  order  of  magnitude. 

In  applying  the  singular  value  decomposition  to  our  problem,  we 
used  the  algorithm  described  in  Reference  11. 

Given  that  we  have  an  inversion  algorithm  we  must  now  consider  the 
following  problem  The  distance  R  between  detectors  in  Eq .  (9)  runs 
over  the  range  (0,°°),  In  order  to  effect  the  inversion  (by  singular 
value  decomposition  or  any  other  method),  one  must  let  the  spacing 

JT  ' 

V.  Blake  and  K.  Barakat,  "The  Inversion  Problem  for  Turn  Feld  Phota- 
r.leatrov  Counting  Statistic (land.  ,7.  rhijfi . ,  M,  19?5, 
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between  detectors  become  infinite,  at  least  in  principle!  Obviously 
this  cannot  be  done  experimentally  so  that  the  question  to  be  answered 
is;  "What  is  the  acceptable  upper  limit  on  b  (call  it  b  1  for  which 
the  inversion  works  in  the  noiseless  case?" 

In  order  to  answer  this  question,  we  first  consider  the  direct  pro¬ 
blem:  given  the  crosswind  profile  determine  1(b).  This  is  simply  the 
integration  of  v(Z)  over  the  path  with  respect  to  the  path-weighting 
function  W(Z,b). 

Before  continuing  these  calculations,  we  choose  represent at ive 
values  of  l.,  1),  and  \ 

1,  =  1.5  •  10  cm  *  1.5  km 


I)  a  1.27  cm 


\  =  6.35 


cm 


Based  on  these  values  we  have 


a  a  0.412 


!•  =  3 . 08  cm 


These  values  will  be  used  for  all  numerical  calculations  in  this  paper. 

The  behavior  of  the  path-weighting  function  W(Z,b)  is  shown  in  lig. 

2  Evaluation  of  W(Z,b)  was  by  high  order  quadrature  sufficient  to 
guarantee  four-digit  accuracy.  As  the  detector  spacing  b  is  increased, 

W  selectively  accentuates  various  parts  of  the  optical  path  with  large 
spacings  concentrating  on  the  end  near  the  laser.  Incidentally,  in  a 
working  system  the  detectors  would  remain  fixed  and  the  apparent  detector 
separation  would  be  done  by  scaling  the  optics. 

IVe  employed 

v(Z)  °  650  e  “  cos(2nZ)  cm/sec  (22) 


(23) 


IS 


*  6S0  e 


ci's  ( 4n Z)  cm/ sec 


as  tost  crosswiml  profiles.  Recall  that  1  knot  SO  cm/soc.  I'ho  resul¬ 
tant  shapes  of  the  slope  of  the  space-time  covariance  function  at  zero 
time  lag  as  a  function  of  detector  separation  b,  L(.b ) ,  are  shown  m  I  i  . 
3.  Needless  to  say  they  display  a  complicated  behavior  for  small  b,  but 
then  tend  slowly  to  zero  in  monotone  fashion  for  large  b.  Large  number 
of  sign  changes  in  v(Z)  manifest  themselves  in  a  more  complicated  be¬ 
havior  of  l;.(b)  for  small  6.  Several  other  crosswind  profiles  were  run: 
they  all  indicate  that  ultimately  1(b)  tends  to  zero  as  b  is  increased. 

Based  on  this  information,  we  set  b  "3  and  r  *4,  then  in 

ma  x  mu  x 

verted  the  profile  given  by  Lq.  (22)  using  a  uniform  spacing  of  Ab  0.1 

(M  =  31)  and  requiring  N  *  21.  The  inversions  are  shown  in  I  ig.  4  for 

b  =3  (open  circles)  and  b  =4  (solid  circles).  Both  vield  excel- 
max  1  max 

lent  fits  to  the  profile  for  Z  >_  0.4;  however  for  Z  *■  0.4  the  inversion 

points  for  b  =3  oscillate  about  the  true  profiles  whereas  those  for 
‘max  * 

b  =4  are  still  very  good.  It  is  an  artifact  of  singular  value  de¬ 
composition  that  the  inversion  points  for  Z  =  0,  1  are  always  zero. 

There  are,  of  course,  21  singular  values  since  N  =  21.  Values  of  b  >4 

*  ma  x 

did  not  yield  reconstruction  significantly  better  than  those  for 

b  =  4.  On  the  other  hand,  profiles  inverted  for  b  '  3  were  almost 

max  1  max 

useless  due  to  wild  oscillations.  At  least  with  respect  to  the  profiles 

we  inverted,  it  appears  that  8m  3  4  is  a  reasonable  compromise  between 

accuracy  and  physical  realizability. 

liven  in  the  "noiseless"  case  we  do  not  employ  all  21  singular  values 
in  the  reconstruction;  rather  we  employ  Id  singular  values  because  '  ,  ^ 

and  o  are  enough  in  error  to  cause  serious  distortions  in  the  recon¬ 
struction  of  the  crosswind  profile. 

In  order  to  mimic  the  experimental  situation,  we  add  signal  depen¬ 
dent  noise  to  U(b)  in  the  following  fashion: 


li(b)  =  (1  ♦  Su)E(b)  .  . 

noisy  no i sol  os « 


(24) 


Here  <5  is  a  small  number  and  u  is  a  random  variable  governed  by  a  roe 
tangular  probability  density  function 


f(u)  =  \  . 


1  <  n  <  1 


0  ,  elsewhere 


(23) 
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100  v(z)cm/sec 


Z 


Figure  4.  Inversion  of  v(z)  corresponding  to  Eq.  (22)  by  singular  value 
decomposition  using  19  o's  in  the  noiseless  situation:  solid 
is  actual  profile,  open  circles  correspond  to  p  =3  and 
solid  circles  to  bmw  *  4.  mdx 

Ilia  X 
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A  typical  inversion  realization  of  the  profile  given  by  liq.  t J J ) 
is  shown  in  Fig.  5.  The  amount  of  "noise"  added  to  li(p)  via  F.q.  (.24) 
was  51  (i.e.,  5  =  .025).  According  to  t he  recipe  quoted  in  lq.  (21), 
this  means  that  we  should  use  lb  singular  values  in  the  reconstruction 
of  the  profile. 

A  second  profile,  Lq.  (23),  was  also  reconstructed  and  the  results 
are  shown  in  Fig.  (».  The  solid  circles  are  the  reconstructed  values  of 
the  profile  in  the  noiseless  case  with  ti  =  4.  The  solid  triangles  are 
the  reconstructed  values  of  a  typical  sample  realization  with  31  noise 
(i.e.,  3  =  .015)  added  to  E(8).  The  results  speak  for  themselves. 

As  we  noted  in  the  introduction,  Heneghan  and  Ishimaru  used  an  in¬ 
version  scheme  dependent  on  statistical  regularization.  It  is  difficult 
to  compare  their  results  since  most  of  their  calculations  are  not  ad¬ 
dressed  to  the  same  problem  with  which  we  are  concerned,  namely  the 
point-wise  inversion  of  v(Z).  However,  in  Fig.  b  of  their  paper  they 

show  an  inverted  velocity  profile  using  data  obtained  from  Harp1'  which 
is  similar  to  that  of  Fig.  5  of  the  present  paper. 


IV.  COMMENTS 

The  previous  calculations  have  demonstrated  that  it  is  possible  to 
obtain  the  crosswind  profile  using  the  correlation  slope  method  via 

singular  value  decomposition,  when  C“(z)  is  constant  over  the  path. 

However,  as  a  reviewer  has  pointed  out,  typical  spatial  and  temporal 

y 

fluctuations  in  C“  along  the  path  can  overwhelm  the  effect  of  IV  (z,.  1 

by  just  varying  receiver  separation.  It  is  suggested  that  one  might 
overcome  some  of  the  problems  bv  using  higher  order  receiver  spatial 

filtering  techniques  such  as  described  in  Clifford,  et  a  l1'1  and  Ochs, 


‘%V.r.  Harp,  "A  If  Pr  :\  'n  Frp, .'i,v 

Motii'iet  -md  Turbulrnt.  Struct ><»v  •  >/  thr  .Ifwv.' •. 
Hadlr8/"ir>u'r  Lib ,  Hnii'.,  PFI - .1  .* .  V'.'. 

14 

S.F.  Clifford,  /. H.  Oi'hs  .nni  W:vd 

'hr  SrrintilLi'-ionx  rf  ■?  h \ndrr-  S, •<?»,',”  Apr'.  •'  ' « 
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Figure  6.  Inversion  of  v (Z )  corresponding  to  Eq.  (23).  Solid  line  is 
actual  profile,  solid  circles  correspond  to  noiseless  situa 
tion  with  8  =  4,  solid  triangles  correspond  to  3%  in 

E(8)  »ith  6^,  *  4. 
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et  al,^  Lee.*b.  Such  an  analysis  has  been  initiated  by  one  of  us  (KB). 
See  also  Leuenberger,  et  al.*' 

It  is  tempting,  when  C“(z)  varies  over  the  path,  to  invert  the  in- 

,  n  ? 

tegral  equation  relating  to  C“  to  the  log-amplitude  of  the  covariance 
for  zero  time  delay  using  the  inversion  scheme  (e.g.,  Reference  18) 

t 

which  will  guarantee  C“  >_  0.  The  resultant  C“  could  then  be  substitute  d 
into  Eq.  (1)  for  the  crosswind  profile.  Unfortunately,  this  approach 
will  fail  because,  as  Strohbehn  has  shown,  one  cannot  distinguish  be- 
tween  a  C“  layer  and  a  change  in  the  turbulence  spectral  power  law. 
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